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Abstract 

A capability for modeling ice crystals and mixed phase icing has been added to GlennICE. 
Modifications have been made to the particle trajectory algorithm and energy balance to model this 
behavior. This capability has been added as part of a larger effort to model ice crystal ingestion in aircraft 
engines. Comparisons have been made to four mixed phase ice accretions performed in the Cox icing 
tunnel in order to calibrate an ice erosion model. A sample ice ingestion case was performed using the 
Energy Efficient Engine (E 3 ) model in order to illustrate current capabilities. Engine performance 
characteristics were supplied using the Numerical Propulsion System Simulation (NPSS) model for this 
test case. 


Introduction 

Mason (Ref. 1) describes a situation where an aircraft engine can encounter rollbacks and flameouts 
at high altitude conditions due to ice crystal ingestion. Numerous in-fight encounters had been observed. 

It was hypothesized that the cause of the incidents was the ingestion of a high volume of ice crystals into 
the engine. More so, this report theorized that the ice crystals could be tolerated if the crystals remained 
small solid particles. However, this paper suggested that the ice crystals first melt upon contact with warm 
engine parts. Subsequent ice particles are then more likely to stick to the wetted surface, allowing ice 
accumulation to occur. The engine abnormalities then occur when this ice buildup sheds. Venkataramani 
(Ref. 2), Veillard (Ref. 3), and Rios (Ref. 4) describe previous attempts to model the engine icing 
problem. This report will describe an overview of the GlennICE and NPSS software. It will then describe 
the modifications made to GlennICE in order to model ice crystal scenarios. Included in these 
modifications is an empirical relationship for ice particle induced mass loss (ice erosion) that used to 
govern mixed-phase conditions. Data for this relationship was taken from a report by Al-Khalil (Ref. 5). 
Finally, a sample case was performed for ice crystal ingestion using the results from the NPSS model to 
provide inputs to GlennICE. 


Nomenclature 


Q drag coefficient 

C p heat capacity 

D ai diffusivity of ice in air 

d particle size 

F force 

H enthalpy 

IWC ice water content 

h heat transfer coefficient 
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k 

thermal conductivity 

L 

characteristic length 

^SLlh 

heat of sublimation 

LWC 

liquid water content 

M 

Molecular weight 

m 

mass flux 

Nu 

Nusselt number 

P 

pressure 

Pv 

vapor pressure 

< 7 

heat flux 

r 

coefficient of restitution 

Re 

Reynolds number 

Sc 

Schmidt number 

Sh 

Sherwood number 

St m 

mass Stanton number 


momentum Stanton number 

T 

temperature 

t 

time 

V 

velocity 

v c 

critical velocity 

P 

collection efficiency 

s 

strain rate 

X 

mass transfer parameter 

fl 

viscosity 

il 

freezing fraction 

P 

density 

a 

shear stress 

Xm 

dimensionless mass parameter 

^mom 

dimensionless momentum parameter 

X/ 

dimensionless thermal parameter 

CO 

mass fraction 

Subscripts: 

a 

air 

conv 

convection 

ev 

evaporation 

f 

freezing 

film 

film 

heat 

heat added 

i 

ice 

KE 

kinetic energy 

m 

melting 

0 

original 

rb 

runback 

s 

surface 

sens 

sensible 

sub 

sublimation 

t 

threshold 

w 

water 

00 

ambient 
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GlennICE Overview 


The GlennICE (Ref. 6) project was initiated in 1999 as a means to develop a consistent set of 
standards and practices for creation, use and maintenance of NASA administered icing computational 
tools. One goal of this process is to provide a mechanism for a more team-based approach to software 
development. Current development of icing codes at NASA has been more separated and individual- 
based. The plan has been to establish a common computational framework for integration of NASA icing 
analytical tools, and to create a set of requirements both for the development of these tools and for their 
integration into the common framework. The features of the resulting system would include providing a 
common user interface for all tools and allowing user selection of analysis modules. A modular design 
would allow more streamlined substitution of new or improved software elements, and the defined 
process would include documented standards for rapid integration of new software. This process also 
requires a version control system designed to allow both control of approved software and flexibility for 
rapid prototyping, as well as documentation that allows rapid identification of software design and use. 

NPSS Overview 

The Numerical Propulsion System Simulation (NPSS) code (Refs. 7 to 9) is used to model the 
thermodynamic cycle of complete aircraft engines and therefore can be effective at predicting the engine 
system operation. The development of this tool is a joint effort between NASA Glenn Research Center, 
the Department of Defense, and Unites States aeropropulsion industry. It is used to model the physical 
interactions between components of a turbine engine reducing the need for full-scale engine tests and 
component experiments ultimately impacting the time it takes to bring an aircraft engine to market. 

At its core the NPSS is an object-oriented code that brings together multiple disciplines such as 
aerodynamics, structures, and heat transfer which supports multi-fidelity (numerical zooming) from zero- 
dimensional to three-dimensional engine component codes. Zooming allows a designer to vary the level 
of detail of analysis throughout the engine based upon the physical processes being studied. It is 
important to consider the aircraft engine as a system of components interacting with each other not simply 
in isolation. In the NPSS these component models, potentially at various levels of fidelity, are allowed to 
interact with each other to see the effects they have on the overall engine performance. A simulation of a 
complete aircraft engine will provide sufficient detail of the individual components and their interactions 
with other engine components over a range of flight operating conditions. This allows for an examination 
of many design options before costly tests are conducted. 

GlennICE Modifications 

Three modifications have been made to GlennICE for handling ice particle impacts. First, there are 
the modifications to the particle trajectory routine that track the energy transfer to and from a particle. 
Second, there is a model for mass loss due to erosion in a mixed phase environment. Finally, there are 
changes to the mass and energy balance whereby the incoming ice particles can add to the ice mass at the 
surface. In addition to these modifications, an analysis was performed to assess if coupling between the 
air and ice phases is necessary. 

There are two separate modifications to the particle trajectory equations in order to handle ice 
crystals. First, a separate drag equation is needed. The drag equation for liquid droplets includes an 
increase in drag due to particle deformation whereas ice crystals remain spherical. It is also assumes that 
ice crystals do not breakup due to shear. In addition, drag correlations were added for ice particles that are 
cylindrical as well as disk-shaped. The second modification was to track the change in particle 
temperature (and phase change if necessary) through the air stream. Previously, GlennICE did not track 
droplet temperature as they were assumed to be in equilibrium with the air. However, the mixed phase 
capabilities have been added to analyze engine ice events. In this scenario, the air, water drops and ice 
crystals will increase in temperature when they enter the engine. Since the air will respond more rapidly 
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to this increase, there will exist a driving force for energy transfer between the particles and the 
surrounding air stream. 


Particle Drag Equations 

The following equations for particle drag were obtained by curve-fitting drag data for spheres, 
cylinders, and disks. The reference used was Perry’s Chemical Engineers’ Handbook (Ref. 10), but 
several references contain the same data. 

For spheres 


For cylinders 


For disks 


24 

C d = 1- 0. 4 + 

d Re 

24 

C d = — + 0.3 + 
d Re 


6 

1 + VRe 

6 

1 + VRe 


for Re < 1 00 
for Re > 100 


C d = 4. 1 94 Re -0 ' 93 1 for Re < 0.0 1 

C d =13.344 Re“ a691 for0.01<Re<20 

C d =3.2344-0.3341og(Re) for 20 < Re < 1 000 

C d = 0.93 for Re >1000 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

(6) 


24 6 

C d = — + 0.4 + = for Re < 40 (7) 

Re 1 + VR^ 

C d = 1.5*10“ 6 Re 2 +0.001Re + 1.8176 for 40 < Re < 1 000 (8) 

Q =1.3 for Re >1000 (9) 

Particle Energy Balance 

The following analysis will present the energy equation for a spherical ice particle. The particle 
energy equation is simply the increase or decrease in enthalpy of the particle balanced with the losses due 
to convection and sublimation or 


dH 

~dt~ 


^conv ‘Zsub/evap 


( 10 ) 


Prior to or after any phase change by melting of the ice crystal, this equation can be given as 
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P i^~ p i , . hud (T a T s ) + nd p a h m L sub (co a co^) 

6 dt 

(ID 

for spheres 

After a particle has started to change phase, L su b is replaced by L em p in the equation above. It will be 
more convenient to solve this equation in non-dimensional form, hence 

-* t hd 

Nu = — 
k a 

(12) 

Sh = h ' nd 

D a,i 

(13) 

p r _ ^P,a^a 

K 

(14) 

Sc = M a 
P a D a,i 

(15) 

M a P 

(16) 

Cp'jPjd~ 

lT ~ 12 k a 

(17) 

The resulting energy equation then becomes 


dT ’ = Nu (t -t )+ Sh Pri “ b (m„-<o s ) 
dt 2 x t v a SJ 2 x t Sc C p i v s ' 

(18) 

where 


For spherical particles 


Nu = 2 + 0.6 Re 1/2 Pr 1/3 

(19) 

Sh-2 + 0.6Re 1/2 Sc 1/3 

(20) 

where 

D P a\ V a~ V i\ d 

Re = — 1 — 

(21) 

For cylindrical particles 
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For discs 


Nu = 0.3 + 0.62 Re 1/2 Pr 1/3 


1 + 


' Re ^°’ 625 


, 0.8 


.2.82 * 10 5 ) 


f ' -> 2/3 3 ^ 


1 + 


2 0.4 2 
v Pr y 


Sh = 0.3 + 0.62 Re 1/z Sc 


l/2o„l/3 


1 + 


' Re ^ 


x0.8 


V 


2.82 * 10 3 


f 


1 + 


''(X^ 

V Sc y 


2/3' 


1/4 


Nu = 0.664 Re 1/z Pr 


1/2 p r 1/3 


Nu = 


0.037 Re 0-8 Pr 


1 + 2.443 Re _ai (pr 2/3 -l) 
Sh = 0.664 Re 1/2 Sc 1/3 


Sh = 


0.037 Re a8 Sc 


l + 2.443Re“ ai (sc 2/3 -l) 

Reduction in particle size for a sphere is given by the mass balance 


(Re <200,000) 

(Re >200,000) 

(Re <200,000) 

(Re >200,000) 


Pi 


^ Kci 3 ^ 


V 6 J 


= Shin/ 2 p a D C(>/ 


d 


or after integration 


d 2 = i / 2 - Xt 


X = 


4Shp a Z) fljf 


(®J CIqq ) 


( 22 ) 


(23) 


(24) 

(25) 

(26) 
(27) 


(28) 


(29) 

(30) 


by 


During phase change, the ice/water temperature remains fixed and the rate of melting of ice is given 


pjLf ^ hud (T (l T m )+nd h m L eyap ((o a to TO ) 


(31) 
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By fixing the crystal temperature at the melting point, the melting rate can be calculated directly from 
this equation. 


Coupling Between Particles and Air 

In previous icing analysis, the airflow and particle dynamics were calculated separately, which 
implies that coupling of the equations are unnecessary. In this regime, it should be clarified whether or 
not coupling is needed. Crowe (Ref. 1 1) estimated mass coupling to be important if 

Z >1 
st M 

(32) 

where 


z = IWC 
P« 

(33) 

cu _ ^ m 

(34) 

Pi d 2 

» yyi / \ 

6 Shp a D a ^ (COj, — COag ) 

(35) 

If we try to maximize the coupling, we need the maximum ice water content, the smallest 
size, the lowest ice velocity, the longest compressor length and the highest mass transfer rate. 
9 g/m 3 , d= 15 pm, F iee = 40 m/s, P = 27,000 Pa, Sh = 5, L = 1.4 m and (eo s -ec>oo) = 0.02 then p, 
z 

0.044 and = 0.018 , thus mass coupling is unimportant. 

^ m 

Momentum coupling is important when 

particle 
Let IWC = 
j 0.4, x m 

z >i 

1 4 - qt 

± -r - J t n om 

(36) 

where 


T V- 

Cf - ' 

^ L mom 

(37) 

T _P dd 2 
V 

(38) 

thus using the previous values, 


x v = 0.000744 

(39) 

St mom =0-021 

(40) 

0 0225 

Z = 0.0225 * = 0.022 

1 + 0.021 

(41) 
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so momentum coupling can be ignored. 
Energy coupling is given by 


Z Lsuh =0.1 (42) 

Cf r T 
Jl m ^ p,i 

using the previous values so energy coupling is unimportant as well. Based on this analysis, the one-way 
coupling performed in icing software such as Glennlce is still justified. 

Coefficient of Restitution 

When a solid particle hits a surface, the resulting collision is considered elastic or inelastic. If the 
collision is elastic, none of the incoming energy is absorbed by the surface. For inelastic collisions, some 
of the energy is absorbed by the surface and some is used to determine the velocity of the particle after the 
collision. This energy ratio is called the coefficient of restitution. If this coefficient is equal to one, the 
collision is considered elastic. If this coefficient is equal to zero, all of the energy is absorbed by the 
surface. In this case, all of the mass must remain on the surface as well. Therefore, in the ice crystal 
environment being modeled, the coefficient of restitution can be used to determine the fraction of 
impinging ice crystals that remain on the surface, similar to the splashing model used for supercooled 
large droplets. Higa et al. (Ref. 12) measured the coefficient of restitution for an ice sphere with a 3 cm 
diameter at velocities up to 0.7 m/s and a wide range of temperatures. They curve-fit the restitution 
coefficient to the following equation: 


r — 




Jcj 



(43) 


They reported values of the critical velocity (V c ) after which particle collisions were no longer quasi- 
elastic at velocities of 25, 35, and 61 cm/s at 269, 261, and 245 K, respectively. The decrease in 
restitution coefficient was caused by fragmentation of the ice spheres. Based on this equation, most ice 
particle collisions at the surface will have essentially a zero restitution coefficient, meaning that almost all 
of the energy of the collision is transferred to the surface and the ice crystals are expected to fragment. 
Initial build up of an ice layer is then caused in part by the kinetic energy of the collision being used to 
melt the surface ice fragments. GlennICE assumes that the fraction of ice crystals present on the surface 
will be (1-r) as this is the fraction of energy transferred to the surface during the collision with the 
surface. However, this does not necessarily mean that the ice crystals will accumulate. In order for 
accretion to occur, a film thickness must be present. Otherwise, the ice fragments would be swept away 
by the flow field. In the limited mixed phase experiments presented below, the ice crystals do not 
accumulate even in a glaze ice condition where a water film is present. However in this case the water 
film is limited to a very small region near the stagnation zone. In the current model, the energy and mass 
balances are calculated assuming the ice crystals remain on the surface. For control volumes that do not 
have a liquid film, the ice crystal mass is removed prior to the geometry modification routines. As more 
experimental data becomes available, this assumption can be revisited. Erosion effects, as described 
below, are accounted for in all regimes. 


Ice Erosion 

Ice erosion occurs when ice particles hit an ice shape that has already accreted. The force of impact 
can dislodge portions of the ice shape, especially individual ice feathers. Al-Khalil (Ref. 5) performed 
tests in the Cox icing tunnel in order to document the mixed phase capabilities of that facility. The model 
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was a 36 in. chord NACA0012 airfoil. All cases were run for 10 min. at 100 kts, a = 0° and a water MVD 
= 20 pm. Results from Al-Khalil’s report show that ice impacts can cause erosion of an ice shape in a 
mixed phase test. Four cases were selected for calibration with the GlennlCE erosion model. The cases 
are given in Table 1. 

The cases at the lower temperature (cases 3 and 4) were unaffected by the ice crystals as the ice 
thickness compares well with values from GlennlCE without any erosion effects. This suggests that the 
low air speed of these conditions is not high enough to cause erosion. At the warmer condition (cases 14 
and 15), the introduction of an IWC of 0.7 g/nT caused a 20 percent reduction in the ice shape as shown 
in Figure 1. The amount of ice loss was calculated using the utility program THICK. This result seems to 
indicate that the lower bonding strength of the ice at the warmer condition is enough to cause the onset of 
ice erosion. Rist (Ref. 13) measured the deformation and fracture of 1.5 cm ice particles over a wide 
range of temperatures. He measured the strain rate at failure to the corresponding stress applied and 
produced the following relationship: 


= (cyi -cr 3 )‘ 


, 4.2 


exp 


14- 


Q_ 

RT j 


(44) 


In this equation, Q is the activation energy that was measured as 69 kJ/mol in the experiment and R is 
the gas constant. For the cases listed, the warmer condition would require only half as much force to 
produce the same strain rate in the ice as the colder temperature cases. The impacting force of each 
particle is proportional to the mass of the particle (and thus proportional to the diameter cubed), the 
velocity and the particle frequency. Without performing a stress analysis of each ice shape there is 
insufficient data to build a proper correlation due to the sparseness of the test matrix. However if it is 
assumed that the cases without erosion effect are at the erosion threshold, the erosion effect can be 
estimated. First, the impact stress was assumed proportional to the force of the ice particle impacts. Thi s 
force was determined by multiplying the particle momentum with the frequency: 


K ,3 


V 

( IWC^ 

— p id 

V 



|_6 


d 

V Pi V 


(45) 


The term in the second bracket is the particle frequency. This equation would make the stress 
proportional to d 2 , V 2 , and IWC 13 . This would make the applied stress the same for the cases where 1WC 
= 0.7 g/m 3 . However, since the temperature is higher in case 15 than case 3, twice as much strain rate is 
produced in the ice, causing a 20 percent reduction in ice size. The model in GlennlCE first assumes that 
case 4 from the Al-Khalil data is at the threshold where erosion becomes important. For this condition and 
for conditions with lower shear stresses due to particle impact, no erosion takes place. For conditions with 
higher impact energy, the ratio of strain rate can be estimated by 


s 

G 


f d ^ 

2 

|V 

2 

' IWC ' 

U'J 


u. 


liwcj 


1/3 


exp| 

f 

exp 

V V 


f u-R- 

v RT 


14-- 


RT, 


tJJ 


(46) 


The erosion threshold values are d t = 150 pm, V, = 100 kts, IWC, = 0.7 g/m 3 , and T, = 267 K which is 
the average ice temperature in the impingement region for this condition. If it is further assumed that the 
erosion rate is proportional to the increase in strain rate, the erosion rate can then be related back to the 
incoming conditions. Comparing case 4 with case 15 from Al-Khalil’s data gives: 
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f 200) 

l 

{ 100^ 

1 

f 0.7 

V 150 ) 


1.100/ 


1.0.7 J 


1/3 


exp 


14 — 


69 

S.3 14 * 1 0 — 3 *273 


v exp l 


14 — 


69 

3.314* 10 -3 * 267 J 


= 3.5 


(47) 


Given that case 4 does not appear to have any erosion effects while case 15 shows a 20 percent loss in 
mass due to erosion leads to the following empirical equation: 


erosion = 0.08 


f — -l 

j 


(48) 


where s t is the strain rate at the erosion threshold. Obviously, much more data would be necessary to 
determine a more accurate value for the erosion threshold as well as to verify the dependence of erosion 
on the icing parameters. 


Surface Energy Balance 

The mass and energy balance add terms for the added mass of ice, the kinetic impact of the ice 
particles, and the melting of ice crystals. For a glaze surface, 

W hmp,w ^irnpj — m f + m rb + m ev + m sub + w fllm (49) 

where m imp ; = (l - erosion )(l -r)p,IWCF (50) 

r is the restitution coefficient calculated earlier. 

The energy equation becomes 

*717111 *7heat *?KE,a *?KE, w *7KE,i — ^conv — *7evap — *7sub — *7sens — 0 (51) 

where qy m is the heat removed by additional freezing or heat added by melting of ice. As described below, 
GlennlCE assumes that only one of these processes can take place in a particular control volume. 

The first term that is added due to ice crystal impacts is the kinetic energy transfer of the impacting 
crystals, 


_ 1 . 

*7KE,i - ~ m imp, i 



(52) 


where r is the restitution coefficient calculated earlier. The second term added is the melting term given 
above. The final term that is changed is the sensible heat transfer from the drop or particle. Previously, the 
drops were assumed to be at the ambient air temperature. Currently, the particle temperature is calculated 
as shown earlier and the resulting temperature at impact is transferred to the icing balance. In addition, 
any phase change or changes in particle and droplet size are accounted for in these equations. 

In the solution procedure, the terms other than melting and freezing are summed. Based on the sign of 
that term, the remaining energy is used to either freeze additional surface water or used to melt incoming 
ice. In a situation where more energy is carried away allowing additional freezing, then the freezing 
fraction of the incoming ice crystals is equal to one and the additional freezing fraction of the incoming 
water is calculated from the energy balance. In a situation where more energy is available to melt the 
incoming ice, then the freezing fraction of the incoming water is zero and the fraction of incoming ice that 
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is melted is again calculated from the energy equation. In any situation where the sum of the other terms 
is exactly zero, then the freeing fraction of the incoming water is zero and the freezing fraction of the 
incoming ice is equal to one. The solution procedure determines the mass of ice or water directly and does 
not rely on calculating a freezing fraction. 


Results 

The energy efficient engine, or E 3 , was designed at NASA Glenn for the intention to lay an advanced 
technology foundation for a new generation of turbofan engines. It was a cooperative program with 
industry aimed at developing and demonstrating advanced component and systems technologies for 
engines. Since the geometry is public domain, it was selected for initial tests of the computer model. 
Figure 2 shows a top view of this engine while Figure 3 shows the cross-section. 

The NPSS code was used to fly the E 3 engine over a typical flight envelope to determine locations in 
the engine where ice accretion might occur. The flight envelope provided the NPSS code with the altitude 
and throttle settings for the engine performance analysis. The NPSS code was used in a 0-D 
(thermodynamic cycle analysis) mode to determine locations in the compression system where ice 
accretion might be supported. Figure 4 illustrates the total temperature variation over a portion of the 
flight envelope for the E 3 at engine station 2.2. Based on icing event statistics of commuter and large 
transport aircraft reported by Mason (Ref. 1), the E 3 was flown over the flight envelope at various altitude 
and temperatures (at and above standard atmospheric conditions). The representative flight envelope 
shown in Figure 4 includes two cruise altitudes and one descent profile. During the two cruise portions of 
the flight, the engine operates at reduced thrust after attaining the desired altitude. At the end of the 
second cruising altitude, the engine is spooled back to 5 percent maximum thrust in preparation for 
descent. The aircraft then descends from maximum cruising altitude to prepare for landing at a reduced 
thrust. An engine operating point at 22 R above standard atmospheric conditions within the 39000 ft 
cruise portion of the flight envelope was selected for testing GlennlCE. Figure 5 shows the fan tip 
performance map of the engine at the elevated temperature. The fan operates well away from the surge 
line over the flight profile. However as indicated in Figure 6, the fan hub and quarter stage compression 
system operates near the surge line and with geometry changes due to ice accretion has the potential to 
push the compression system toward surge thus compromising the engine performance and causing a 
safety concern. 

The initial analysis performed for this report focuses on the region from the inlet through the low- 
pressure compressor (FPC) to the entrance of the high-pressure compressor (HPC). The system modeled 
in NPSS using a 0-D analysis was used to obtain temperature, pressure, and flow rates at various engine 
stations. A diagram of the region analyzed is shown in Figure 7 and the results of the NPSS analysis at the 
selected location in the flight envelope are shown in Table 2. 

A zero-order analysis only provides values at specific locations such as the conditions before and 
after the fan for example. It does not include the dimensions of the system, hence the use of the term zero- 
order. While NPSS is capable of performing a higher fidelity model, the zero-order results were sufficient 
for the purposed of testing the initial model. In order to use this information in GlennlCE, a grid was 
created in GRIDGEN and the values from Table 2 were used to supply velocities, pressures, and 
temperatures at the grid points. Since a flow analysis was not performed using this grid, it can be 
somewhat coarser than that used for Navier-Stokes analysis since its only purpose is for tracking particle 
trajectories. There were 10 grid blocks used with a total of 19,000 grid points in order to interpolate the 
NPSS results to a grid that GlennlCE could use. Spherical ice crystals of 200 pm MVD were then tracked 
through this artificial flow field to find collection efficiencies on the two splitters. This is shown in 
Figure 8. 

The secondary splitter had minimal impingement and no ice growth, so results for that geometry are 
not shown. An icing time of 5 min and an IWC of 8 g/m 3 were used to generate the results for this case. 
Velocity, temperature and pressure were taken as the conditions at the engine inlet as shown in Table 1. 
Results from the energy balance and mass balance are shown in Figures 9 and 10, respectively. Figure 1 1 
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shows what the ice shape would look like without erosion effects. Based on the preliminary ice erosion 
model used, this ice shape would be completely eroded by the impact of other ice crystals due to the much 
higher velocities inside the engine. It should be emphasized that this result has not been experimentally 
verified. 


Future Work 

This model must be considered at the early stages of development at this juncture. Much work 
remains in order to produce a model ready for release. In particular, experimental data is needed to verify 
erosion effects and the deposition rate of ice crystals when a water film is present. Just as early research in 
the supercooled droplet regime provided references of droplet splashing at low velocities that did not 
coincide with results obtained in the NASA icing tunnel at higher velocities, the models presented in this 
report for coefficient of restitution, particle breakup and ice erosion at low velocities need to be confirmed 
at the higher velocities inside the engine. In addition, while particles impinge on the lower wall of the 
input geometry (shown in Fig. 7), GlennICE needs to be modified to allow deposition and accretion on 
this surface. Experimental data on actual engine geometries in an ice crystal environment is also 
necessary for validation. It would also be useful to test the model with two-dimensional flow solutions 
from NPSS. Once a viable model is ready, it will also be implemented into LEWICE3D so that the full 
geometry can be analyzed. 

Conclusions 

A preliminary model for multiphase ice accretion was developed for GlennICE. An analysis of typical 
ice crystal inputs showed that one-way coupling could still be used in the icing model. The effects of 
temperature change, phase change and sublimation of the ice crystal were added to the trajectory model. 

In addition, non-spherical ice crystal capability was also added. A factor for the coefficient of restitution 
was also added. The energy equation was revised to include the effects of ice crystals, including a 
preliminary model for ice erosion. Erosion results were calibrated to a mixed phase icing test performed 
by Al-Khalil in the Cox icing facility. Icing results were produced for a test case using the E 3 geometry. 
The results show that while ice could form on the splitter lip based on the energy balance, it would be 
quickly eroded by other ice particle impacts. However, this result has not been experimentally validated. 
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TABLE 1.— MIXED PHASE ICING RUNS 


Case no. 

T 

J Ol 

°F 

MVD, ice, 
pm 

IWC, 

g/m 3 

LWC, 

g/m 3 

3 

12 

150 

0.3 

0.7 

4 

12 

150 

0.7 

0.3 

14 

22 

N/A 

0 

0.7 

15 

22 

200 

0.7 

0.7 


TABLE 2.— FLOW CONDITIONS FROM NPSS ZERO-ORDER ANALYSIS 


Station 

T, 

R 

P, 

psia 

Mach 

V, 

mph 

w, 

lbm/s 

1 

462.2 

4.27 

0.78 

776.4 

1313.3 

2 

462.2 

4.27 

0.57 

580.2 

1313.3 

2T 

462.2 

4.27 

0.57 

580.2 

1017.5 

2H 

462.2 

4.27 

0.57 

580.2 

295.8 

21 

513.3 

5.94 

0.48 

520.0 

224.0 

22 

530.7 

6.67 

0.48 

534.6 

202.8 

23 

530.7 

6.67 

0.48 

534.6 

96.6 

24 

530.7 

6.67 

0.21 

240.0 

106.2 

25 

530.7 

6.67 

0.40 

439.8 

106.2 

13 

527.4 

6.47 

0.46 

502.5 

716.3 

14 

527.4 

6.47 

0.48 

526.1 

716.3 

15 

527.4 

6.50 

0.48 

531.4 

813.0 
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Figure 1. — Erosion effect with mixed phase inputs. 



Figure 2. — Top view of energy efficient engine. 
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Figure 3. — Cross-section view of energy efficient engine. 


dT = +22 R 



Altitude, ft 


Figure 4. — Total temperature range at engine station 2.2 of the E 3 over a prescribed 
flight envelope at 22 R above standard atmospheric conditions. Engine operating 
point in cruise at 39000 ft was selected for testing GlennICE. 
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Figure 5. — Fan tip performance map showing corrected flow rate versus total pressure ratio, 
for a range of corrected speed lines. Key engine operating points along the vehicle trajectory 
are superimposed onto the fan map. Corrected speed lines are expressed as a percent of 
the design rotational speed. 
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Figure 6. — Fan hub and quarter stage performance map showing corrected flow rate versus 
total pressure ratio, for a range of corrected speed lines. Key engine operating points along 
the vehicle trajectory are superimposed onto the fan map. Corrected speed lines are 
expressed as a percent of the design rotational speed. 
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Figure 7. — Engine region analyzed with NPSS and GlennICE. 
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Figure 10. — Mass balance terms for E test case. 


NASA/TM— 2011-216978 


18 



NASA/TM— 2011-216978 


19 


REPORT DOCUMENTATION PAGE 

Form Approved 
OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 
data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 
burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 
Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 
control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 
01-02-2011 

2. REPORT TYPE 

Technical Memorandum 

3. DATES COVERED (From - To) 

4. TITLE AND SUBTITLE 

Mixed Phase Modeling in GlennICE With Application to Engine Icing 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 



5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

Wright, William, B.; Jorgenson, Philip, C.E.; Veres, Joseph, P. 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

WBS 609866.02.07.03.04 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
John H. Glenn Research Center at Lewis Field 
Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

E- 17618 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S 
ACRONYM(S) 

NASA 

11. SPONSORING/MONITORING 
REPORT NUMBER 

NASA/TM-201 1-2 16978 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified-Unlimited 
Subject Category: 03 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

A capability for modeling ice crystals and mixed phase icing has been added to GlennICE. Modifications have been made to the particle 
trajectory algorithm and energy balance to model this behavior. This capability has been added as part of a larger effort to model ice crystal 
ingestion in aircraft engines. Comparisons have been made to four mixed phase ice accretions performed in the Cox icing tunnel in order to 
calibrate an ice erosion model. A sample ice ingestion case was performed using the Energy Efficient Engine (E 3 ) model in order to 
illustrate current capabilities. Engine performance characteristics were supplied using the Numerical Propulsion System Simulation (NPSS) 
model for this test case. 

15. SUBJECT TERMS 

Turbojet engines; Aircraft icing 

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION OF 
ABSTRACT 

uu 

18. NUMBER 
OF 

PAGES 

25 

19a. NAME OF RESPONSIBLE PERSON 

ST1 Help Desk (email:help@sti.nasa.gov) 

a. REPORT b. ABSTRACT c. THIS 

U U page 

u 

19b. TELEPHONE NUMBER (include area code) 

443-757-5802 


Standard Form 298 (Rev. 8-98) 
Prescribed by ANSI Std. Z39-18 



































